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1. Introduction 

Recent more stringent emission regulations have pushed for the development of 
more fuel-efficient and low-emission combustion systems for both ground-based and 
flight gas turbine applications. Experimental estimate of pollutant formation is limited 
due to the inability to probe non-intrusively, the highly three-dimensional flame zone. As 
a result, detailed understanding of where (and why) the pollutants (primarily, NOx, CO 
and unburned hydrocarbons, including soot) are forming is difficult, if not impossible to 
determine using experimental data. This is a major stumbling block in the design of the 
next generation, fuel-efficient gas turbine engines. Numerical prediction could help but 
most codes currently in use employ steady state methods and are unable to capture the 
unsteady dynamics that control the pollutant formation process. In addition to the 
inaccuracy in the prediction of NOx formation, CO emission is also not predicted very 
well under dynamic conditions. For example, recent studies (Bhargava, 2000) have 
shown that as the equivalence ratio is reduced, CO emission first decreases and then 
suddenly increases rapidly. This sudden increase occurs during lean combustion when 
perturbation to the flame can cause instability and also lead to incomplete combustion. 
The condition at which this sudden increase occurs is very important to predict and 
control since it impacts the performance and efficiency of the gas turbine engine. 
However, so far, no numerical simulation tool has been able to predict this behavior. 

Accurate prediction of unsteady combustion requires a comprehensive model that 
can predict flame structure and propagation characteristics, pollutant formation and 
transport, and ignition/extinction phenomena over a wide range of flow conditions in 
high Reynolds (Re) number, three-dimensional (3D) flows as in a gas turbine combustor. 
Direct Numerical Simulations (DNS) are not practical since the resolution and 
computational resource requirements far exceed the present and even future capability. A 
simulation approach that has become popular in recent years is large-eddy simulation 
(LES). In LES, all scales larger than the grid resolution are numerically simulated using a 
space- and time-accurate scheme while the effect of scales below the grid resolution is 
modeled using a subgrid model. Although many LES studies have been reported, a 
validated LES approach for practical systems is yet to be demonstrated although the 
approaches discussed in this paper have shown a serious potential to achieve this 
capability. 
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Over the last few years, there has been a major effort at NASA-GRC to develop a 
unique computational capability called the Numerical Propulsion Simulation System 
(NPSS). In this approach, a simulation model that will allow simultaneous calculation of 
the flow in the entire engine (i.e., inlet, compressor, combustor, turbine and nozzle) is 
currently being pursued and developed. Within this structure, the simulation tool for the 
combustor, called the National Combustion Code (NCC) has reached a major state of 
maturity and is currently being evaluated for practical applications both at NASA and at 
industrial R&D sites (e.g., GEAEC). Various advanced models have been implemented 
to make NCC a state-of-the-art modeling tool and its application to various problems are 
being demonstrated and/or discussed in numerous recent papers (e.g., Liu, 2001, Raju, 
2001, Lannetti et al., 2001, Moder, 2001, Ajmani, 2001, Shih et al., 2001). 

The NCC is based on Reynolds-Averaged Navier-Stokes (RANS) modeling 
approach and is the current choice since it is computationally efficient when complex 
systems have to be simulated. However, it has some inherent (and well-known) 
limitations. In RANS modeling, the effect of all turbulent length and time scales are 
modeled using global turbulence models such as the popular k-e model. Closure of the 
turbulent reaction rate and turbulent species flux are also based on global model that 
ignore some fundamental physics of turbulent mixing and combustion. For example, the 
dynamical interactions between flow, unsteady heat release, acoustic (pressure) 
oscillation and turbulence is lost since only the time-averaged flow field can be simulated 
with such an approach. These very important phenomena directly impact fuel-air mixing 
process and are controlled by processes that occur at the small-scales. Unfortunately, 
even unsteady RANS models cannot resolve these fine-scale interactions in a physically 
consistent manner since all these effects are modeled. 

In this three-year research study, we plan to take advantage of two crucial 
developments: (a) the demonstrated capability of NCC to simulate flows in complex 
configurations using robust algorithms and (b) the recent demonstration of a unique 
subgrid combustion model for LES applications. We (the team from Georgia Tech and 
NASA/GRC) will combine these two capabilities into a new simulation model, denoted 
hereafter NCCLES that will contain all the capabilities of the original NCC but will now 
contain additional algorithms that will allow it to be used in a LES mode. It is envisioned 
that NCCLES will eventually become a specialized tool that will address issues that the 
original NCC cannot. 

2. Technical Objectives 

The proposed research is carrying out both model development and model 
implementation simultaneously. The model development effort (which will primarily be 
at Georgia Tech) is addressing issues related to complex high-Re flows, including 
supersonic flows. Model implementation of the subgrid kinetic energy model and the 
subgrid mixing and combustion model within NCC is being carried out at NASA/GRC. 
Substantial progress in both the collaborative effort and in the new advances in the LES 
subgrid models has been achieved already. More details of the progress made to date are 
summarized below. 
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The overall technical objectives of this three-year project are: 

1 . Collaborate with the NCC team and provide algorithms that will allow them to 
convert the current NCC algorithm to a LES version. 

It is worth noting that since many of the baseline algorithms are already 
developed, transition of this technology to the NCC team has already begun. More 
details of the current collaborative effort are summarized below. 

2. Further develop the subgrid models for momentum and scalar transport for 
implementation into NCCLES. 

This is a major effort since NCCLES is expected to deal with many 
different scenarios and the subgrid models have to be able to simulate all such 
flows without compromising the physics. Therefore, the models currently 
operational need to be re-evaluated in this context. This effort has already begun 
and we have made substantial progress to-date. 

3. Collaborate with the NCC team to demonstrate and validate the next generation 
NCC for LES applications. 

All the models developed during this research will be implemented in our 
in-house code, LESLIE3D that has been extensively used for such studies. Once 
the validation of the new algorithms is completed, we will transition this 
technology in a timely manner to the NCCLES for further evaluation. Close 
collaboration with the NCC team will be maintained during this effort. 

3. The LESLIE3D Code 

The LESLIE3D code is a fully compressible, 3D finite volume solver that has 
been used extensively in the past to carry out LES studies of complex turbulent reacting 
flows as in the gas turbine combustor. The code can be used for both DNS and LES and 
is second-order accurate in time and has an option to do either second or fourth-order 
accurate scheme in space. A constant coefficient model and a localized dynamic subgrid 
model (where the model coefficients are obtained as a part of the simulation and hence, in 
this case, the closure has no parameters to set or adjust) for the subgrid kinetic energy are 
options in the code. In addition, models for subgrid combustion are also implemented for 
both premixed and non-premixed combustion. For completeness, the LES model 
equations and closures are summarized here. 

3.1 The LES Equations 

The compressible LES equations are obtained by spatial filtering of the 
conservation equations of mass, motion, energy and species. The final set of LES 
equations are: 
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Here, the filtered total specific energy E is given by: E = e + + k sgs , where 

k sgs is the subgrid turbulent kinetic energy given by k sgs = (u,u, - u l u l )/2 . The unclosed 
terms T- gx , H sgs , cr. gs , <f> f pi , © ^ represent the subgrid stress tensor, heat flux, viscous 
work, species mass flux and species diffusive mass flux, respectively, and are: 
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The pressure p is obtained from the equation of state for a thermally perfect gas, 
p = pRT . The mixture gas constant is given by R = ^ Y m R m . A similar expression also 

m 

applies for mixture specific heats C p and C v where C p -C V =R. 

3.2 Subgrid Closure for Convective Fluxes 

Closure of the subgrid flux terms is achieved using a transport equation for the 


subgrid kinetic energy per unit mass k sgs = 




(Menon et al., 1996): 
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The two model coefficients, C v and C £ , are 0.067 and 0.916, respectively for the 

constant coefficient model based on earlier studies (Chakravarthy and Menon, 2001a, 
2001b). A localized dynamic approach developed using a scale-similarity approach can 
also be used to obtain these coefficients model dynamically (Kim and Menon, 1999, Kim 
et al., 1999, Kim and Menon, 2000). Given the subgrid kinetic energy, the unclosed terms 
in LES equations can be modeled. For example, the sub-grid heat flux and the subgrid 
mass flux are: 


H sgs 

O sgs 

i,m 


= ~P 


v, dh 

Pr, dx : 


= ~P 


SC t DX; 


Where Sc t and Pr f are the turbulent Schmidt number and Prandtl number. At present, 
these are assumed to be unity. 


LESLIE3D contains closure models for both premixed and non-premixed mixing 
and combustion process. These are briefly summarized below. 


3.3 Subgrid Flame Speed Model 

This approach is only applicable to premixed combustion in the flamelet regime. 
In this case, the flame is much thinner than the smallest scales of turbulence and 
therefore, turbulence can only wrinkle the flame but cannot affect the flame structure. 
Flame wrinkling and stretching increases the effective flame surface area and therefore, 
can increase the burning velocity. In the flamelet limit, a simplification can be carried out 
whereby the flame is modeled as an infinitely thin surface that separates the unbumt 
reactants from the burned products. In this approach, reaction kinetics need not be solved 
and the problem essentially reduces to the tracking of the flame front. This approach has 
been used quite successfully to simulate lean premixed combustion in full-scale gas 
turbine engine such as the General Electric LM6000. We will be using this version of the 
LESLIE3D initially to develop and evaluate the NCCLES. 

The model involves the transport solution of a progress variable G that lies in the 
range [0,1] with the flame located between the two limits. The G- field is convected by 
the flow field and also propagated normal to itself with the normal burning velocity. In 
the laminar case, this is the classical laminar burning velocity (which contains all the 
thermo-chemical effects of finite-rate kinetics) but in the turbulent field, a closure is 
needed to estimate the effective turbulent burning speed. The LES model for the flame 
front is obtained by the solution of the following filtered G-equation: 
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Here, the terms of the right-hand-side require closure. Here, S sgs is the turbulent 
flame propagation term typically given as: S sgs = S’t’IvgI where S T = S T (u' ,S L ) is the 


turbulent flame speed that is modeled in terms of the laminar flame speed and the subgrid 


turbulence (obtained from subgrid kinetic energy: k sgs = — 



). In the present 


study, a model proposed by Pocheau has been used to obtain the turbulent burning 
velocity. The second term G { sgs ={u i G — u i G) is the subgrid convective flux that also 

requires closure. Various models have been proposed for this term: (a) localized dynamic 
eddy diffusion model (Kim and Menon, 2000), (b) flame curvature model (Peters, 2000), 
etc. We have evaluated all these models as reported below. 

Heat release is incorporated by including the heat of formation in the definition of 

the internal energy using the following relationship: e = C v T +e j H(G - Gq ) where, 


H(Q) is the Heavy-side function and G 0 is an arbitrarily chosen value in the range [0,1] 

where the flame is located. More details are given elsewhere (Kim et al., 1999, Kim and 
Menon, 2000) and therefore, avoided here, for brevity. Although this approach is 
computationally very efficient (since species equations and/or finite-rate kinetics do not 
have to be modeled), this approach cannot be used to capture species variations and/or be 
used to predict pollutant formation. To do this, one has to resort to the more complex 
finite-rate kinetics closure. 


3.4 Subgrid Scalar Closure using Eddy Diffusivity 

For both premixed and non-premixed combustion, if detailed kinetics is of interest 
multi-species equations have to be solved along with the LES equations. In this case, two 

terms, the sub-grid scalar flux, = p[V i m Y m -V i m Y m ] and the filtered reaction rate 
term, d> m require closure. In a conventional approach, a subgrid gradient approximation: 

pv'VY 

@im ~ L — ~ where Sc m is a turbulent Schmidt number, is employed for the subgrid 

Sc m 

scalar flux. Note that, since large-scale motion is resolved in a LES, associated counter- 
gradient processes are resolved (even when a gradient closure is employed for 0-^). The 

closure for <b m s complicated due to its highly non-linear nature. Assumed PDF subgrid 

models, subgrid Eddy Breakup and subgrid Eddy Dissipation models, all of which are 
variants of the models used in RANS codes, have been proposed for this closure. 

However, at this time, the accurate closures for the reaction rate within the 
conventional LES context have yet to be established. In order to ensure accuracy, we 
have developed a new closure model based on a subgrid simulation approach. This is 
described in the next section. 
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3.5 Subgrid Scalar Closure using LEM 


In the LEM closure, the species fields evolve within each LES cell due to 
localized, stochastic reaction-diffusion processes, and are then transported across the LES 
cells due to convective flux. The local sub-grid domain within each LES cell is resolved 
on a one-dimensional domain that resolves all scales of motion and thus, processes within 
this domain can be considered a localized 1-D Direct Numerical Simulation. Thus, all 
scales down to the Kolmogorov and Batchelor's scale are resolved in this simulation 
(Menon et al., 1993, Calhoon and Menon, 1996, Chakravarthy and Menon, 2001a, 
2001b). 

To describe this method, consider an “exact” reaction-diffusion equation for a 

scalar 


d(j> _ 

— -+ («, + u i ) 
ot 


80 

dx t 



80_ 

dx 


] + ® $ 


( 5 ) 


Here, u i and u\ are the resolved and unresolved velocities, respectively. Also, 
u\ is due to both the Lagrangian convection through the cell face u' lag , and turbulent 
convection at scales smaller than the resolved grid, u' stir 

The one dimensional linear-eddy domain, in every LES cell is initialized with a 
constant number of LEM cells. The number of linear-eddy cells that are needed to resolve 
all the sub-grid scales is computed based on a subgrid Reynolds number probability 
density, obtained over the whole of the computational domain. The length of the linear 
eddy domain is set equal to the LES filter width. Since the Lagrangian convection across 
the LES cell face is based on the mass, the volume of the linear-eddy cell is also carried 
as an additional scalar. At the start of the simulation, volume of each linear eddy cell is 
set to have an equal fraction of the LES cell volume. At each LES step, after updating the 
density, momentum and the energy field, LEM is used to evolve the scalar field. Note 
that, as opposed to the conventional scalar closure approach described in the previous 
section, no scalar equations are solved on the resolved grid. 

Using a fractional splitting technique, we can split the equation (5) as follows: 

Subgrid Processes: f [ u \ tir + -A- [ D . ^-] + cd Adt (6) 

f ox , ox, ox 

t o 1 1 

f -r , d<f> 

Supergrid Convection: ~ — («• + u lag )— ( 7 ) 

ixt les OX; 

Equation (6) represents the sub-grid reaction-diffusion processes that occur 
locally within each LES cell using the LEM model. Three processes occur within each 
LES cell: (a) turbulent convection at scales smaller than the resolved grid (which is 
modeled using stochastic stirring events called triplet maps, (b) molecular diffusion and 
finite-rate kinetics, and (c) volumetric heat release. Additional details are given in cited 


NASA/CR— 2003-2 12202 


7 



references. After the first two subgrid processes, volumetric gas expansion caused by 
heat release is implemented by expanding each linear-eddy cell in the domain by an 

n 

amount equal to AV^ = . 

Pi 

Equation (7), which represents the advection (convective transport) of the scalar 
field by the resolved velocity field is modeled by a volume of fluid (VOF) approach 
(Chakravarthy and Menon, 2001a, 2001b) in which the linear eddy cells are convected 

Ar 

across LES cell face at the advective time scale given by A t conv = (CFL)mm — — . 

Ideally, splicing in a compressible simulation can be performed every A t conv step. But, in 

order to provide a close coupling between the fluid dynamic and the flame front 
evolution, splicing is performed every LES time step. This considerably improves the 
numerical stability of the simulation. 

On a general three-dimensional grid the mass transfer from any control volume to 
any neighboring control volume is also predominantly 3D. As a result, different numbers 
of linear-eddy cells are transported in different spatial directions. In the numerical 
implementation of the scalar advection, the 3D advection operator is approximated by a 
sequence of three, one dimensional advection operators. So the order in which these 
operators act on the scalar field can have a significant effect on the scalar field evolution. 
In the new algorithm, the linear-eddy cells going out in the direction of largest out-flux 
are fluxed out first from the right end of the sub-grid domain. Similarly the largest influx 
is added first to the left end of linear eddy domain. This upwind type flux update has been 
used successfully in the past (Chakravarthy and Menon, 2001a, 2001b). 



Arrows 

- Convection Direction 

Color Lines 

- Flame 

Blue Cell 

- LES Cells 

Black Cell 

- LEM Cells 


Figure 1. Schematic illustrating the Splicing Algorithm used for scalar convection. 

3.5.1 New Splicing Algorithm Development under this Project 

The earlier studies were limited to low-Mach number flows and therefore, some 
issues were substantially simplified and easy to implement. For example, in low-speed 
flows, the CFL time step is the convective time step, which is also the same time needed 
for the scalar transport via splicing. However, in the present study, we are eventually 
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interested in extending this approach to fully compressible flows, including supersonic 
flows. In this case, the CFL time-step is the acoustic time step whereas the scalar 
transport is still restricted by the convective time-step, which is at least an order of 
magnitude larger than the acoustic time-step. Thus, there is a possibility of de-coupling 
between the fluid dynamics (which evolves at the acoustic time) and the scalar transport 
(that occurs at the convective time). A substantial effort has been spent on understanding 
this issue and we have now developed a new convective transport model applicable for 
compressible flows. 

The present splicing algorithm advects the sub-grid scalar field based on the 
mass, as opposed to the volume of the LEM cell that was used in the past. The advantage 
of using mass instead of volume is that, it ensures mass conservation at the resolved and 
unresolved scales at each instant of the calculation. Also in the earlier approach, partial 
LEM cells that have to be convected are accumulated, until they reach an integer value. 
In the new algorithm, every time a partial cell has to be convected out, a single cell of the 
same volume and mass, as that of the partial cell and with all the scalar properties is 
convected out. This potentially avoids the mass conservation errors and spurious 
numerical instabilities in the computation. 

Volumetric expansion due to heat release and splicing by itself, also lead to a non- 
uniform distribution of the volume of the linear-eddy cells. To avoid programming 
complexities in a parallel environment, it is desirable to have same number of linear eddy 
cells everywhere in the computational domain. So, to equalize the volume of the linear- 
eddy cells and to retain the same number of linear eddy cells everywhere in the 
computational domain, the linear eddy domain is re-gridded to have cells of equal volume 
after each splicing. Errors introduced due to the spurious diffusion associated with re- 
gridding have not been evaluated yet. But it is expected to be minor, since at most only a 
couple of cells are affected. 

Furthermore, to reduce any biasing error introduced by the operator splitting, the 
order in which reaction, diffusion and stirring are performed is switched after every 
splicing. An additional advantage of this switching is that, any discontinuities introduced 
by splicing are smoothed quickly by turbulent stirring, thus preventing any spurious 
diffusion. 


3.5.2 Subgrid Simulation using LEM 

The details of the sub-grid one-dimensional simulation using LEM are briefly 
summarized here. LEM is a stochastic model that treats reaction-diffusion and turbulent 
convection separately but concurrently. Reaction-diffusion processes evolve on a one- 
dimensional (ID) domain in which all the characteristic length scales in the turbulent 
field (from the integral scale L to the Kolmogorov q) are fully resolved (6 cells are used 
to resolve q). The orientation of the ID domain is in the direction of the scalar gradient 
(and thus, for premixed flame, is in the flame normal direction), and within this domain, 
the equations governing constant pressure and adiabatic flame propagation are: 


dt 


1 dpY t V t 
p dx 


Ml 
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+ F, 


kstir 


( 8 ) 


NAS A/CR— 2003-2 1 2202 


9 



( 9 ) 
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The equation of state for the scalar mixture is P = pT'S' - - and the caloric relation is 

hw k 

n r r 

given by h k = Ah, k + C k (T')dT' . Here, T is the temperature, P is the thermodynamic 
pressure, R u is the universal gas constant and p is the mass density. Also, 
Y k ,W k ,C pk ,d> k ,h k ,V k and A h° f k are respectively, mass fraction, molecular weight, specific 

heat at constant pressure, mass reaction rate, enthalpy, diffusion velocity and standard 
heat of formation (at standard temperature, T°) of the /c-th species. The mixture-averaged 
specific heat at constant pressure and thermal conductivity are respectively, C p and k . 

Fickian diffusion velocity law: V k = , where Dk is the mixture-averaged 

Y k dx 

diffusivity of the /c-th species, is used for molecular diffusion. 

The convective terms u i dY k /dx i and u i dT jdx i in the species and the energy 


equations are represented as F^stir and F Ts tir, respectively. These terms are implemented 
using stochastic re-arrangement events called triplet maps, each of which represents the 
action of a turbulent eddy on the scalar fields. It has been shown that this mapping can 
capture correctly the physical increase in scalar gradient (without affecting the mean 
scalar concentration) due to eddy motion. Three parameters are needed to implement 
these turbulent stirring events: the typical eddy size /, the eddy location within the ID 
domain and the stirring frequency (event rate) X. The eddy size in the range L to q is 
determined randomly from an eddy size distribution, f(l) which is obtained using inertial 
range scaling in three-dimensional turbulence: f(l) = (5/3) l' 8/3 / (r/' 5/3 - L 5 ' 3 ). Here, q is 
determined from inertial range scaling law q = N n L Re' 3/4 where N t) is an empirical 
constant. This constant reduces the effective range of scales between L and q but does not 
change the turbulent diffusivity, as described earlier. 

The event location is randomly chosen from a uniform distribution and the event 
rate (frequency per unit length) is determined by: 


54 i/Re [(L/t7) 5/3 -1] 
TC x L 3 [1 - (/// A) 43 ] 


( 10 ) 


The time interval between events is then given as At st i r = 1/(X,X LE m) where X L em is the 
length of the ID domain. 


The unique feature of the LEM model is that although the scalar evolution is 
simulated in ID, the effect of turbulence on the scalar fields is modeled using 3D scaling 
laws. As a result, flame wrinkling occurs at spatial and temporal scales that mimics the 
effect of realistic 3D turbulent eddies on the flamelets. This formulation has two 
constants: C*. and Nt, both of which arise from the use of scaling laws. These parameters 
are obtained by comparing LEM predictions to experimental data in the flamelet regime 
(Smith and Menon, 1996). The present study uses these same values (C>, = 15 and 
Nr, =4). 


NAS A/CR— 2003-2 1 2202 


10 



4. Progress to Date 


Here, we summarize the progress made to date. It is worth noting that the project 
officially started at Georgia Tech in late January 2002 and therefore, the progress 
noted below is for a period of 9 months. We report here primarily on the effort carried 
out at Georgia Tech over this period. There were two major areas of progress: (a) 
transition of LES technology to NCC team at NASA/GRC and (b) further development of 
LES model for more complex flows. 

4.1 Collaboration with NCC Team at NASA/GRC 

Discussions with the NASA/GRC team members were carried out to set up the 
work plan. It was determined that one area that needs to be explored is the manner by 
which LES is carried out since it differs substantially from RANS application. It was also 
decided that this “learning” process needed to be carried out prior to implementing the 
LES routines into NCC. In order to do this Georgia Tech provided the NCC team with 
our advanced LES code, LESLIE3D (Large-Eddy Simulations using Linear Eddy in 
3D). This code is the primary research tool employed at Georgia Tech to develop all the 
algorithms to be eventually implemented in NCC. Close collaboration has been 
established in order to “teach” the method for doing LES using our code. 

4.1.1 LES Simulation strategy using LESLIE3D 

In order to establish an understanding of how to do LES, the first phase of 
collaboration was to use the LESLIE3D to simulate a test case that has also been 
simulated at Georgia Tech. The NCC researchers have been using the LESLIE3D code to 
carry out the baseline simulation of the GE LM6000. This test case was also simulated at 
Georgia Tech and the detailed analysis package was also provided to NASA/GRC. 

The test case is a full-scale GE LM6000 that was experimentally studied at 
GEAEC. Figure 2 shows the test configuration. 


Sideview 


"Blocks 0.003m thick 



End view 



0.07m 


* 


Figure 2. Schematic of the GE LM6000. 
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Premixed methane-air combustion is simulated under highly swirling conditions. 
The Reynolds number based on the inlet mean velocity and diameter is the 350,000, the 
inlet Swirl Number at the dump plane is 0.56, the inlet temperature is 688 Kelvin, the 
combustor mean pressure is approximately 6.5 ATM. We have employed relatively 
coarse grid of around 600,000 in the past to simulate this problem using the localized 
dynamic subgrid kinetic energy model and the dynamic flame speed model. A typical 
result is given in Fig. 3 which shows the centerline velocity decay. Various closure 

models for the subgrid diffusion term G sgs = (ujG - ufi) were attempted but the results 

are relatively insensitive to them and all cases show reasonable agreement with data. 
Additional comparison of the velocity profiles also show good agreement. 



Figure 3. Centerline mean axial velocity variation. Symbols are the experimental 
data from GEAEC. The various lines correspond to the subgrid closure for the 

convective flux term G sgs = (ii fi - UjG). All models show reasonable agreement 
with data. 


4.1.2 DNS of Decaying Isotropic Turbulence 

Another area where collaboration has been established is in the area of inflow 
turbulence and initialization process. Typically, isotropic turbulence field is added to the 
mean flow. A version of the LESLIE3D was also developed and provided to NASA/GRC 
that can be used to generate isotropic turbulence. Since this solver has the ability to do 
both DNS and LES it was decided to simulate this sort of flow using both the second and 
fourth-order spatial schemes so that a baseline data base can be established to evaluate 
the ability of the NCC base scheme to simulate isotropic turbulence. 

One case chosen for this purpose is direct numerical simulation (DNS) of 
decaying isotropic turbulence. Since the physics of the isotropic flows (such as the 
energy decay rate, form of the spectrum, etc.) are well established from the numerous 
experimental and numerical studies from the past, this exercise can serve as a powerful 
test-bed for establishing the capabilities of the code. 
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Decaying isotropic turbulent flow experiments of Comte-Bellot and Corrsin 
(CBC) (1971) were chosen for validating the compressible solver. In the experiments, 
measurements were taken at three locations downstream of the turbulence-generating 
grid. The mesh size of the grid was 5.08 cm and the mean velocity ahead of the grid was 
10 m/s, giving a mesh Reynolds number of Re = 34000. Downstream of the grid, the 
wind tunnel had a slight contraction (1.27:1) that served to provide an isotropic condition. 
Conceptually, homogenous isotropic turbulence decays in time. However, in the wind 
tunnel environment, the turbulence decays in space as it evolves downstream. In order to 
convert a spatial evolution to the temporal evolution, Taylor’s hypothesis is used and the 
turbulence in space is related to the decay in time of an imaginary box of homogenous 

P dx — 

turbulence through the relation: t = — — , where V (x) is the mean velocity of the 

o U (x) 

flow at the streamwise location x. 

To design a computation for the simulation of the CBC experiment, appropriate 
initial conditions that has the same spectrum as in the experiment is needed. In general 
the initial velocity field for the isotropic turbulence must satisfy at least three conditions: 
(a) conservation of mass, (b) the generated velocity field must be a real function of space 
and time, and (c) the turbulent field must have a realistic energy spectrum. 

The first two conditions are required in order to obtain a stable numerical solution 
while the third conditions is needed to reduce the initial transient as the velocity field 
evolves to form realistic isotropic turbulence. In reality, not only does one need to 
specify the energy spectrum, but also the relative phases of the modes. In practice, 
however, this phase information is not available so that fields are usually constructed 
with random phase, but in such a way that they have a prescribed initial energy spectrum. 
Therefore, this initial condition will lead to a transient in which the phases adjust 
themselves to appropriate values. 

The governing conservation equations for compressible flows are solved using a 
finite volume scheme that is second order in time and second or fourth order in space. In 
general, for a given Reynolds number, the second-order scheme requires twice the 
resolution of the fourth-order scheme to maintain DNS and/or LES accuracy. This is an 
important issue for model and code validation. The results reported below are for DNS 
and LES using the fourth-order scheme. The computational grid used is 64 x 64 x 64 with 
uniform spacing (for DNS) and periodic conditions are imposed on all the boundaries. 
The parameters corresponding to the simulation and the experiments are given in Table 1. 


Table 1 


Station 

Ut 

L 

Re* 

Simulation 

Grid 

1 

42 

71.6 

LES 

32 x 32 x 32 

2 

98 

65.3 

DNS 

64 x 64 x 64 

3 

171 

60.3 

DNS 

64 x 64 x 64 
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The DNS simulations are initialized with a velocity field corresponding to the 

Re A =65.3 (station 2, — =98) and is evolved for a time corresponding to the equivalent 

L 

spatial distance between station 2 and station 3. Figure 4 shows the energy spectrum 
along with the experimental spectrum at two of the locations. It can be observed the 
present solver captures the energy cascade in the inertial range very well. The most 
significant aspect of this figure is that dissipation scales are resolved accurately. This is 
very important for any solver to demonstrate. If the numerical dissipation of the scheme 
is too high, dissipation range cannot be resolved very well. On the other hand, if the 
dissipation is not sufficient, the code becomes numerically unstable. The fact that the 
Kolmogorov cascade and the dissipation range of the spectrum are captured accurately is 
a proof that the code provides the correct amount of dissipation, while still retaining the 
stability of the solution. The deviation in the largest scales is because of the smaller size 
of the physical box used for the simulation. Since a bigger box size would require a 
higher resolution, slightly smaller box is used in the current simulation. But, as it is 
pointed out earlier, it is the smaller scales that are usually difficult to capture and hence, 
this should not be of any concern. 

Figure 5 shows the vortical structures at the initial time and at a time when the 
turbulence has fully evolved. The decay of the turbulence and the growth of the small 
scales is clearly apparent in this figure. 

4.1.3 LES of Decaying Isotropic Turbulence: 

Having established the accuracy of the numerical scheme using the DNS, it 
remains to establish the accuracy of the sub-grid models employed in the present 
simulation. For this purpose, a large eddy simulation with a higher Taylor’s Reynolds 
number ( Re^ = 71.6) was carried out. The initial velocity field was initialized with a 

velocity field that would provide the correct resolved and subgrid energy (as that in the 
experiments) corresponding to the computational grid (32 x 32 x 32) employed and the 
Re^. A one-equation sub-grid kinetic energy along with an eddy- viscosity type closure is 

used to model the sub-grid scale processes. The rationale behind this and the validity of 
this approach has been very well established and documented. 

Figure 6 shows the decay of the turbulent energy decay rate. As can be seen the 
current model accurately captures the energy decay rate from the large scales to the small 
scales. Sub-grid models are designed to provide dissipation from the resolved to the sub- 
grid scales and it is apparent that the current subgrid model provides the right amount of 
dissipation since otherwise it would not be possible to capture the correct decay rate 

Figure 7 shows the energy spectrum from the present LES and the CBC 
experiment. Inertial range is captured very well despite the coarse grid used in the 
simulation. As expected, the energy spectrum deviates from the experiments at the small 
scales, since the dissipation range scales are not resolved explicitly in any LES. Finally, 
Figure 8 shows the evolution of the vorticity field from the current LES at the initial time 
and at a later time. 
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Figure 4. Energy Spectrum for Re^ = 65.3 corresponding to CBC experiment. 



Figure 5. Vorticity Field corresponding to Re A = 65.3 at t = 0 and t = 0.2 seconds. 
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Figure 6. Kinetic Energy Decay in time. 



Figure 7. Turbulent Kinetic Energy spectrum for Re A = 71.6. 



Figure 8. Vorticity Field in LES for Re A =71.6 at initial time and at a later time. 
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4.2 Research Development at Georgia Tech 

In addition to the collaboration with the NCC team reported in Section 4.1, 
research was undertaken at Georgia Tech to further develop the LESLIE3D code. Four 
areas have been under focus. These are (a) Flame-turbulence test problem for NCCLES 
validation, (b) LES of supersonic spatially evolving shear layers, (c) LES of combustion 
in the Trapped Vortex Combustor (TVC), and (d) DNS and LES of two-phase mixing 
layers. The reasons for the choice of these test cases are noted in the following sections. 

4.2.1 Flame-Turbulence Test Problem for NCCLES Validation 


The interaction of a premixed flame with turbulence is a well-established test 
problem used in many past studies to understand the impact of turbulence on flame 
dynamics. This problem is particularly suited for code validation since both DNS and 
LES have been carried out in the past (including detailed kinetics) and therefore, the 
physics of this problem is well established. This test case is also relatively simple to 
implement and test. This test problem can also be used for evaluating the algorithm’s 
accuracy in reacting flows. Some results are reported below on the current effort. 

DNS and LES of the turbulent premixed stoichiometric methane-air flames with 
finite-rate effects have been performed. To reduce the cost of computations incurred due 
to the high resolution DNS and by the inclusion of finite rate chemistry, the simulations 
are restricted to two-dimensions. However, this is not a problem in this study since 
similar studies have been performed in the past (Echkeki and Chen, 1991) to understand 
flame-turbulence interactions. 

The governing conservation equations are solved using a finite-volume scheme 
that is second-order accurate in time and fourth-order in space. The DNS was initialized 
in a two-dimensional computational domain with a plane laminar premixed flame. The 
reactant composition is a stoichiometric methane-air mixture at atmospheric pressure and 
temperature. The initialization for the scalars (species & Temperature) is obtained from a 
one-dimensional steady unstrained computations using PREMIX. A single step global 
chemistry (Westbrook and Dryer, 1971) with 5 species is used. The boundary conditions 
are periodic in the transverse direction and non-reflecting at the inflow and outflow 
(Poinsot and Lele, 1991). A computational grid of 400 x 400 with uniform spacing is 
employed for DNS and 100 x 100 for LES. 

The initial turbulent field is generated from a specified energy spectrum, 


satisfying continuity. The spectrum is of the form: E(k) = 16m' -exp-J - 2 — 


where k is the wave number, k 0 (= 6) is wave number of the most energy containing 
scale and u' a is turbulent intensity which is equal to 10.57 m/s in the present simulation. 
The turbulence velocity field is super-imposed on the laminar flame. The ratio of the 
turbulence intensity u[ to the laminar flame speed, S L is 4.2 and the ratio of the integral 

scale to the laminar flame thickness, is 3.4. The turbulence Reynolds’ number based on 
the integral scale is 137 and the Taylor’s micro-scale is 36. 
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Three sets of simulations were carried out to contrast and compare the 
performance of the models employed in the present computations. DNS, LES-LEM and 
another LES using Eddy-Break-up closure (EBU) for scalars were performed. LES-LEM 
and LES-EBU simulations typically employed a uniformly spaced 100 X 100 
computational grid. LES-LEM used 18 LEM cells per LES cell. 

Figure 9 shows the instantaneous temperature contours computed using the LES- 
EBU, DNS and LES-LEM simulations, respectively. The flame topology indicates the 
presence of highly wrinkled flamelet due to flame-turbulence interactions. Note that LES- 
EBU diffuses the flame zone considerably, whereas LES-LEM actually makes the flame 
sharper when compared to the DNS. This is not surprising since LES-LEM actually 
captures the flame within the subgrid with very high resolution (actually better than 
DNS). For example, where as the DNS is capturing the flame over approximately 15 grid 
points, LES-LEM is capturing the flame over approximately 3-4 LES grid points. 
However, since there are 18 LEM cells per LES cell, the actual resolution is over 
50 subgrid points. Thus, LES-LEM captures the flame very accurately. This is similar to 
what was observed earlier (Chakravarthy and Menon, 2001a, 2001b). This is 
demonstrated in Fig. 10 where the instantaneous vorticity and density contours for the 
three cases simulated are compared. (Contours not shown at the same time. Note, this is 
only a qualitative comparison). Vigorous interaction between the local vorticity and the 
flame are apparent from this figure. 



(a) Eddy-Break-up (b) DNS (c) LES - LEM 

Figure 9. Instantaneous temperature contours during flame-turbulence interaction. DNS is 
using 400x400 where as the LES is using 100x100. LES-LEM captures the flame as a sharp 
front (even better than DNS) since the subgrid resolution is better than in the conventional 
DNS. LES-EBU on the other hand diffuses the flame and fails to resolve the wrinkles. 
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(a) Eddy-Break-up (b) DNS (c) LES - LEM 

Figure 10. Instantaneous density and vorticity contours during flame-turbulence 
interactions. Comparison is not shown at the same time. LES-LEM approach captures the 
flame very sharply compared to LES-EBU, even though same number of grid points are 
used in both the approaches, 


4.2.2 LES of Supersonic Mixing Layers 

One of the eventual goals of this project is to extend the NCCLES capability to 
supersonic combustion applications. However, in order to achieve this capability, the 
LES-LEM approach has to be extended to highly compressible flows. 

The goal of the present study in the first year is to extend LES formulation to the 
highly compressible flows and to investigate the effects of the compressibility on the 
flows of engineering interest. In the first stage of this study the dynamic subgrid scale 
models are formulated with explicit correction terms, which accounts for compressibility 
effects. Simulations are underway for spatially evolving supersonic mixing layer for 
which there is some experimental data for model validation. Also, this particular flow 
was investigated earlier (Nelson and Menon, 1998) using an earlier version of the 
LESLIE3D code. 


LES Equations for Compressible Turbulence: 

The governing equations are obtained from Navier-Stokes equations by using 
standard filtering instead of Favre filtering as done in the previous studies. While it 
results in somewhat more complicated form of the governing equations than using Favre 
filtering, it provides for explicit terms that can be identified as “compressible” terms. 
Therefore, the modeling of these terms can be considered independently of the modeling 
issues typically used for the closure of the subgrid shear stresses. 
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The filtered LES equations can be written as 


dp dpu t d 

-di = -^--d^ (PUi - pUi) 


(11a) 
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The resolved viscous stress tensor can be written as t„ = 2 u(—( — - H ) o„ 

lJ 2 y dxj dx/ 3 dx k " 

The total energy is E = e +—u.u. +k sgs , where subgrid kinetic energy is defined as 

2 

1 


k sss = (u i u i - </, </, ) . The filtered equation of state p = pRT + R(pT - pT) completes 
2 

the set of the governing equations. 

Compressibility correction terms: 


The above LES equations contain additional terms that require closure. The 
density-velocity correlation terms in the continuity equation is purely a compressible 
term and this term is rewritten as the difference between Favre and conventionally 

filtered velocity: c- gs = pu t - pu t = p(u - u) . This term would vanish in the in viscid 

limit and it is generally negligible in the region of mild gradients (Chen et al., 1989). In 
the region of stronger gradients, it can be argued that the contribution of this term would 
be proportional to the mean density gradient, which lead to a gradient diffusion model for 

rs 

this term: c- gs = pu t - pu t « —v c — — where v c is a “compressible viscosity” that must be 

dx t 

determined. 

Since this term is expected to be significant only near strong density gradients, it 
is advantageous to formulate it in such a way that it would include a density gradient 
switch, avoiding therefore the excessive dissipation in the smooth mean flow regions. 

From dimensional consideration this term can be formulated as: v c = a c S p \ u k n k | A as 
“compressible viscosity correction term,” where 
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Here, Re A = 


I "k n k I A 

v + v, 


, where A the local grid spacing is the characteristic length scale 


and velocity normal to cell face is the velocity scale. The expression for the scaling 
coefficient a was obtained from the numerical simulations of 1-D non-linear Burger’s 
equation (Nelson and Menon, 1998). 

The turbulent eddy viscosity, used in the cell Reynolds number definition is 

defined as v , « c v 4lP gs A (following Schumann, 1976) and used for modeling of 

‘incompressible’ terms of subgrid scale closure equations. 

The LES continuity equation can now be rewritten as 

dp 

dt 

The subgrid term in the momentum equation can be rewritten as T-f s{c) = pu^u } - pu i u j . 

Again, as in continuity equation, compressibility effects are significant only near the 
regions of high gradients, and again the compressibility effects model should be triggered 
by the presence of high gradients (e.g., shocks). Since shocks were observed to have a 
dissipative effect on turbulence, the compressibility effect is modeled as an additional 

dissipation term M.f s(c) » -2 pv c (S t] - — S kk 8 - ) where resolved rate of strain tensor is 
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) . Thus, the modeled momentum equation can then be written as 
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The modeled subgrid kinetic energy equation in the compressible case is written as: 
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(15) 

This subgrid model contains two model coefficients which can be dynamically 
computed as before (Nelson and Menon, 1998). Finally, the compressibility effects in the 
energy subgrid model closure equation can be estimated as 


dpE 

dt 
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This equation introduces yet another subgrid scale modeling coefficient, which 
can also be dynamically modeled as well. 

Supersonic Spatial Mixing Layers: For validation of the numerical scheme and 
compressibility correction terms the comparison with Samimy et al (1991) experimental 
data was conducted. The original experiments were conducted in a high Reynolds 
number aerodynamic tunnel. Two stream mixing layers at different Mach numbers was 
separated at the entrance to the test region by a splitter plate. The flow scheme is 
provided in Figure 11. The optical test section IS 500x80 mm and velocity measurement 
were conducted using two-component laser Doppler velocimetry. 



Figure 11. Spatially evolving supersonic mixing layer setup of Samimy et al. (1991). 


The incoming flow parameters are provided in the Table 1. 


T (K) 

P, (kPa) 

M, 

m 2 

M c 

u. 

(m/s) 

U/U, 

p2/pl 

5 (mm) 

Re 

291.0 

314.0 

1.80 

0.51 

0.52 

479.5 

0.355 

0.638 

8.0 

9.4x1 0 6 


The computational area starts immediately after the edge of the splitter plate and 
has a length of 400mm and height 152 mm. Both upper and lower walls are slip walls. 
The periodic boundary is used in the z-direction. The inflow is fully specified on the 
supersonic part of the flow while on the subsonic part a characteristic boundary condition 
is used. The outflow boundary conditions are calculated using characteristic form of the 
governing equations using Poinsot and Lele method (1991). The initial conditions for the 
incoming boundary layers are computed using a 2-D boundary layer solver based on the 
known experimental data. For better comparison with the experiment a pseudo-turbulent 
velocity fluctuation are added to the inflow boundary. The turbulence with the assumed 
spectrum is calculated in a box and convected into the inflow boundary, scaled according 
to boundary layer solution on inflow boundary. 

Figures 12 and 13 illustrates the instantaneous density, pressure and span wise vorticity 
contours iso-surfaces for the developed mixing layer obtained in a recent simulation 
using the current LESLIE3D with the compressible correction. A grid 120x91x32 is 
employed for this simulation. This study is still underway and more results will be 
reported in the near future. Figure 14 shows the results obtained in an earlier study 
(Nelson and Menon, 1998), which compares the typical results for the shear growth rate 
and the normalized velocity profiles. 
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Figure 12. Iso-surfaces of density and pressure for the supersonic 

shear layer mixing. 



Figure 13. Spanwise vorticity in the supersonic shear layer. 




Figure 14. Comparison of the mixing layer growth and mean velocity profile with 
the experimental results. These results were obtained earlier (Nelson and Menon, 
1998) but are being revisited for LESLIE3D validation in the supersonic regime. 


NASA/CR— 2003-2 12202 


23 



The next step would be to attempt to compare the numerical simulation results 
with the DNS data for mixing layers with significant compressibility effects. This will be 
an area of research over the next few months. 


4.2.4 DNS of Two-Phase Flows in Temporal Mixing Layers 

Another area of research was in the development of DNS/LES methods for 
simulating single and two-phase mixing process in a temporal mixing layer. This test case 
is also chosen for NCCLES validation for some very important reasons: 

(1) The current NCC code has the capability to do periodic conditions in only 
one-direction. Thus, with this capability, the isotropic turbulence test case 
noted in Section 3 will be quite difficult to accomplish unless major revision 
of the code is carried out. On the other hand, the temporal mixing layer can be 
simulated using periodic in one-direction and slip boundary conditions in the 
other two directions. 

(2) The temporal mixing layer has been used extensively in the past for both DNS 
and LES studies and there is a large amount of data available for sort of flow. 
Also, temporal mixing layer is a subset mixing flow field of direct relevance 
to reacting flows. 

(3) The gas phase temporal mixing layer simulation capability was developed 
many years ago at Georgia Tech. Here, we have extended its capability to 
simulate two-phase flows (e.g., spray mixing). This is with the focus for the 
eventual application to spray combustion, which is the major area of interest 
in the application of NCCLES. 

In the following figures we show some samples of the DNS and LES studies using the 
temporal mixing layers. Results of the present DNS have been compared with DNS 
predictions in the past using a spectral code (Ling et al, 1998). Very good agreement with 
past studies is observed in all cases. The effect of particle size (as determined by the 
Stokes number) was the principal focus these studies. It can be seen that the simulations 
are able to capture all these effects accurately. The effect of using second-order and 
fourth-order spatially accurate schemes using the same grid resolution is also studied 
here. It can be seen that both schemes are reasonably accurate; however, the fourth-order 
scheme is better. In the present case, the second-order scheme is quite good since these 
simulations primarily focused on the rollup and pairing of large-scale coherent structures. 
These coherent structures are primarily 2D in their nature (and as initialized here) and 
therefore, the second-order scheme is reasonably accurate. However, if the mixing layer 
is initialized with 3D isotropic turbulence then it would become apparent that the fourth- 
order scheme would be superior for the same resolution. 
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Figure 15. Evolution of the temporal mixing layer (shown here is the vorticity) 
during vortex pairing in the presence of particles for Stokes number of 4. 



Figure 16. Particle motion and entrainment during vortex pairing (for a Stokes 
Number of 4). Top sequence is from DNS of Ling et al (1998) while the bottom 
sequence is from the present DNS study. 
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Figure 17. Comparison of past DNS (Ling et al, 1998 top series) prediction of 
particle motion with present DNS (bottom series) results. Figures from left to right 
are for Stokes Numbers of 1, 4 and 100, respectively at the same non-dimensional 
time. 




Figure 18. Modal energy and momentum thickness growth rates in the temporal 
mixing layer. The effect of second and fourth-order accuracy of the DNS is shown. 



Figure 19. RMS of the particle number density in a cell as a function of Stokes 
Number is shown. Comparison with past DNS studies is shown. 
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5 . 


Conclusions and Future Plans 


The results reported here show that substantial progress has been made in the 
present research in this project. In addition to providing LESLIE3D codes and helping 
NASA/GRC researchers to understand LES methodology and its application, we have 
been developing additional features that will eventually find application in the NCCLES. 
In particular, three areas of development and validation are underway: (a) Flame- 
Turbulence Interactions using premixed combustion, (b) Spatially evolving supersonic 
mixing layers, and (c) Temporal single and two-phase mixing layers. The configurations 
chosen are such that they can be implemented in NCCLES and used to evaluate the 
ability of the new code. 

The studies planned in the near future for the first year’s effort is to complete the 
LES-LEM application in the flame-turbulence interactions and to further carry out LES 
of the supersonic spatially evolving shear layers. In addition, the extension of the LES- 
LEM approach for spray combustion will be addressed. In the second year, we anticipate 
further extension of the LES approach to more complex supersonic flows including shock 
waves. 

In addition, in the second year we will begin LES study of the Trapped- Vortex 
Combustor (TVC). This combustor was chosen as a final test case for NCCLES and 
General Electric Aircraft Engine Company team had agreed to collaborate with the 
NASA/Tech team in this research. Therefore, we will begin discussion with GEAEC to 
identify a generic TVC configuration of their interest and begin setting up the LESLIE3D 
for this configuration. It is anticipated that some new insight into the TVC flow field will 
be achieved during the second years’ effort under this project. 
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